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disk scale height, and unperturbed surface density in order to study the long- 
term accretion and final masses of gas giant planets. We first derive an analytical 
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I— 1 1 We construct an analytic model for the rate of gas accretion onto a planet 

^ embedded in a protoplanetary disk as a function of planetary mass, disk viscosity, 

O ■ 

I formula for surface density profile near the planetary orbit from considerations 

I of the balance of force and the dynamical stability. Using it in the empirical 

^ I formula linking surface density with gas accretion rate that is derived based on 

! hydrodynamic simulations of Tanigawa & Watanabe (2002, ApJ 586, 506), we 

^ I then simulate the mass evolution of gas giant planets in viscously-evolving disks. 

I We finally determine the final mass as a function of semi-major axis of the planet. 

^ I We find that the disk can be divided into three regions characterized by different 

I processes by which the final mass is determined. In the inner region, the planet 

I grows quickly and forms a deep gap to suppress the growth by itself before disk 

I dissipation. The final mass of the planet in this region is found to increase with 

^ ' the semi-major axis in a similar way to the mass given by the viscous condition 

I for gap opening, but the former is larger by a factor of approximately 10 than 

the latter. In the intermediate region, viscous diffusion of the disk gas limits the 
gas accretion before the planet form a deep gap. The final mass can be up to the 
disk mass, when disk viscous evolution occurs faster than disk evaporation. In 
the outer region, planets capture only tiny amounts of gas within the lifetime of 
the disk to form Neptune-like planets. We derive analytic formulae for the final 
masses in the different regions and the locations of the boundaries, which are 
helpful to gain a systematic understanding of the masses of gas giant planets. 

Subject headings: methods: analytical — methods: numerical — solar system: 
formation 
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Introduction 



A fundamental but unresolved issue with planet formation is how the mass of a giant 
planet is fixed. In the solar system there are four giant planets, which are characterized 
by their massive hydrogen/helium envelopes. Jupiter and Saturn consist mostly of hydro- 
gen/helium, while Uranus and Neptune are mostly of ice but have significant amounts of hy- 
drogen/helium. The giant planets are different in mass: Jupiter's mass is 1 x 10^^ Mq{= Mj), 
Saturn's mass is ~ 0.3Mj, and Uranus' and Neptune's masses are ~ O.OSMj. Furthermore, 
the extrasolar planets detected so far also range in mass from about one Neptune mass to 
~ 10 Mj0; those planets are believed to have massive hydrogen/helium envelopes like the 
giant planets in the solar system. 

There are two competing ideas for the formation of giant planets, the core accretion 



model (e.g..lMizuno 



198d: lBodenheimer &: Pollacklll986l ) and the disk instability model (e.g. 



Camerorull978l : lBosslll989l ). In the core accretion model, a rocky/icy core first forms through 



collisional aggregation of planetesimals, followed by envelope formation due to substantial 
accretion of gas from the circumstellar (protoplanetary) disk. In the disk instability model, 
density fluctuation of the disk gas grows to form a gaseous planet, followed by core formation 
due to sedimentation of heavy elements in its interior. The advant ages and disadvantages 
of both models are discussed in several literatures (e.g., iBosd |2002| ). which is not repeated 
here. In this paper we consider giant planet formation in the context of the core accretion 
model. 

In the core accretion model, the process of the accumulation of the envelopes is divided 
into two phases, the subcritical- accretion and supercritical- accretion phases, in terms of the 
dominant energy source. The transition from the former to the latter occurs when the mass 
of a core reaches a critical value. In the subcritical-accretion phase, incoming planetesimals 
supply energy to the envelope so that the envelope is in the hydrostatic equilibri um. The core 
accret ion thus controls the gas accretion in this phase. The phase 2 found by iPollack et al. 
( 119961 ) is a part of the subcritical-accretion phase. In the supercritical-accretion phase, 
the energy supplied by planetesimals is insufficient to keep the hydrostatic structure of the 
envelope, so that the envelope substantially contracts and releases its gravitational energy, 
resulting in runaway accretion of the disk gas. 

The supercritical-accretion phase can be further divided into two subphases. In the 
former subphase, the gas accretion is con trolled by contraction of the envelope and thus oc- 



curs on the Kelvin- Helmholtz timescale (IBodenheimer fc PoUackl 1 19861 : iPoUack et al.l Il996 



^http:/ /exoplanet.eu/ 
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Ikoma. Nakazawa. &: Emoril l2000l : llkoma fc Gendall2006l ). However, because the gas accre- 
tion driven by the envelope contraction accelerates rapidly with time, the supply of the disk 
gas inevitably becomes unable to keep up with the demand of the cont racting envelope. Thus, 
in th e latter subphase, the disk-gas supply limits the gas accretion (ITanigawa &: Watanabe 
2002! . hereafter TW02). 



In this paper we focus on the gas accretion in the latter subphase of the supercritical- 
accretion phase, i.e., the phase in which the Kelvin-Helmholtz contraction of the envelope is 
under way but the gas accretion is limited by disk-gas supply. After the onset of the runaway 
gas accretion (i.e., in the supercritical-accretion phase), the planet always experiences the 
limited gas supply (i.e., the latter subphase) by the end of its growth. Although there is an 
idea that the masses of Uranus and Neptune were fixed in the subcritical-accretion phase 
because of dissipation of the disk gas (jPoUack et al.lll996l ). we will explore another possibihty 
that the masses of those planets were fixed in the supercritical-accretion phase. 

One idea is that the mass of a giant planet is completely fixed when the planet opens 
a gap in the disk. From a physical point of view, the gap opens when the gravitational 
scattering by the planet overwhelms both viscous diffus ion (the viscous conditio n) and push 



due to pressure gradient (the thermal condition) (e.g.. iLin fc Papaloizoulll993l ). In normal 



proto planetary disks, the final mass of a giant planet is determined by the thermal condition 



[e.g.. Ilda fc Linll2004l ). which is equivalent to the condition that the Hill radius is equal to 
the disk scale height at the planet's location. The final mass is thus determined only by 
temperature of the disk gas for given stellar ma ss and planet 's location. In protoplanetary 
disks similar to the minimum- mass solar nebula (jHayashilll98ll ). the thermal condition yields 
the final mass of approximately one Jupiter mass at 5 AU, but predicts that the masses of 
giant planets increase with the semi-major axis, which is in clear contradiction to the current 
configuration of the solar system. Furthermore, because, in reality, the gap is not vacuum 
but a low-density regio n, continuous gas accretion t hrough the gap takes place, as suggested 



by several studies (e.g.. lArtymowicz &: Lubowlll996l ). The effects of the subsequent accretion 



on the final masses of giant planets have been poorly understood. 

Another idea is that gas accretion is truncated by depletion of the disk gas before gap 
opening (e.g., viscous dissipation and photo-evaporation of the disk gas). This means that 
the final mass of a giant planet is the total mass of the disk gas that the planet has captured 
until the disk gas dissipates completely. This may be responsible for the current configuration 
of the solar system, since it takes longer time for outer giant planets to accrete principally 
because of their long orbital periods. However, there have been no quantitative studies of 
this possibility. 



The purpose of this paper is to gain a systematic understanding of how the final mass of 



-4- 



a giant planet is determined. Specifically, we clarify (1) how much the mass of an accreting 
giant planet increases via slow, continuous accretion through the gap, (2) when and where 
disk dissipation dominates gap opening, and finally (3) how the final masses of giant planets 
depend on several disk parameters. To do so, we simulate the long-term accretion of a 
giant planet in a self-consistent fashion. The essential ingredients of the modeling are global 
viscous evolution of the disk and the flow pattern of the accreting gas, as described below; 
especially, the latter effect was completely neglected by previous studies. Furthermore, we 
derive analytical formulas for the final masses of giant planets and the boundaries between 
several different regimes. 

During the viscous evolution of the disk, gas radially migrates from outer regions to 
the site of giant planet formation. The mass flux limits gas accretion onto giant plan- 
ets in some situations. Because of limited CPU speed of computers, azimuthally-averaged 
one-dimensional simulations of the disk evolution are inevitable to follow the long-term ac- 
cretion. Several previous studies simulated the radial mass fl ux in viscously ev o lving disks 
and then calculated g as accretion rate onto giant planets. iLecar fc Sasselovl (120031 ) and 

Guillot fc Huesd (l2006l ) also simply assum ed that gas accretion rate ont o giant planets is al- 

ways p roportional to the viscous mass flux. IVeras fc Armitagd (120041 ) and lAlibert. Mordasini. fc Benz 
(120041 ) used an approximate formula for gas accretion rate that fits the result s of two- 



dimen siona l hydrodynamic simulations of disk gas with an embedded planet done by lLubow. Seibert. &: Art 
(119991 ) and iD'Angelo. Henning. fc Kleyi (120021 ). However, the formula was derived based on 
a limited number of numerical simulations. Hence, it is in question whether the formula 
applies in a variety of situations. 

To understand the detailed pattern o f accretion flow and to obtain the gas accretion 
rate onto the planet, several workers fe.g..lKlev 



D'Aneelo. Henning. k KlevI 12001 iBate et al. 



2003 



19991: iLubow. Seibert. fc Artvmowicz 



D'Angelo. Kley. fc Henning] l2003l ) per 



1999 



formed two- or three-dimensional hydrodynamic simulations of interaction between a planet 
and disk gas. Since they simulated global disks, the effect of gap formation on accretion rate 
was automatically included. However, because of the time-consuming simulations, they did 
not simulate the growth of giant planets on timescales of the viscous disk evolution. 

This study is clearly different from the previous studies in that this study includes the 
effects of the flow pattern of the accreting gas in simulating long-term growth of planets 
embedded in globally evolving disks. To do so, we use the empirical relation between gas 
accretion rate and unperturbed surface density that was derived by TW02 (see §3 for the 
details). They carried out detailed two-dimensional hydrodynamic simulations of gas accre- 
tion flow onto a planet. They demonstrated the importance of the local flow pattern on the 
accretion rate. 
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In this paper, we first derive an analytical formula for radial surface density distribution 
of protoplanetary disks with an embedded planet in §2. Next, we describe how to include 
the effect of the flow pattern of accreting gas, following TW02, to calculate gas accretion 
rate onto giant planets in §3. In §4, we simulate the long-term accretion of giant planets, 
based on the prescription given in §§2 and 3, and obtain the final masses of giant planets 
for wide ranges of several disk parameters. Also, we derive analytical approximate formulas 
for the final masses of giant planets and the boundaries between several regimes. Finally, we 
conclude this paper in §5. 



We consider a planet embedded in a gas disk, both of which are assumed to rotate with 
the Keplerian velocity around a star. We derive an approximate formula for azimuthally- 
averaged surface density as a function of radial distance to the central star, r. We adopt a 
local coordinate system in which all quantities except the surface density are independent of 
r. We assume steady states and no mass flow toward the central star. We do not consider 
migration of the planet. The validity of these assumptions are discussed in §4.41 

We employ a local coordinate system co-rotating with the planet on a circular Keplerian 
orbit around the central star (a local shearing-sheet approximation). The origin of the 
coordinate system is the planet's position. The x axis is on the line from the star through 
the planet; the direction of the y axis is the same as that of the velocity vector of the planet. 
The equation of motion of the disk gas in the y direction in this coordinate system is given 



where S is surface density, and Vy are velocities in the x and y directions, u is viscosity 
coefficient, and Vy is force acting on gas per unit mass. Here we have omitted the advection 
term in the y direction because of axisymmetry. We also eliminate all the terms on the 
left-hand side of this equation on the assumptions of steady states [d/dt = 0) and no mass 
flow {vx = 0). Hereafter we consider only x > because of symmetry. 

To derive an analytical formula for we assume that Vy is equal to the Keplerian 

velocity, Vy = —{3/2)QpX where Qp is the Keplerian angular velocity at the planet's position. 
This assumption will be found to be reasonable later. Equation ([1]) is thus transformed into 



2. Surface density 



by 






(91nS 

^ dx 



'V - 



0. 



(2) 



The gravity of the planet perturbs the motion of the disk gas. We follow the impulse 
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prescription (e.g., iLin &: Papaloizoulll979l ) that approximates scattering of the disk gas as 
small-angle scattering of a particle by a point-mass object. On the prescription, Vy is given 
by 



97r VmJ "^p^^p"^ 



(3) 



where is the distance of the planet from the central star and Mp and M^, are the masses 
of the planet and the central star, respectively. This approximation is inappropriate for x < 
min(2rH, r^^J{Mp/ M^) /h)^ where th = rp(Mp/3M^)-'^/^ namely, the Hill radius, because this 
region is what is called the horse-shoe region, in which gas is not pushed away from the planet. 
However, since the gas in the horse-shoe region has no contribution to the accretion f low o f 
interest (see in we neglect such a region. Note that ICrida. Morbidelli. &: Massed (120061 ) 
proposed a more elaborate formula for the torque (a ^/-direction force in our formulation) 
that includes the radial transfer of torque through waves, which they found in their two- 
dimensional hydrodynamic simulations. However, since the formula is unable to apply for 
wide ranges of the parameters, as they mentioned, we do not use their formula in this study. 



Substituting equation ([3]) into equation ([2]), we obtain 

'X\ ^3 



S(x) = Soocxp 



Svis(2^)) 



(4) 



where Sqo is the surface density at infinity (i.e., far away from the planet's orbital radius) 
and £ is defined as 



V 



8l7r V 



0.146 



10-5r2fip 



-1/3 



1 1/3 



Mp 
10-3M, 



2/3 



P- 



(5) 
(6) 



This surface density profile is basically the same as that derived by iLubow Sz D'Angelo 
( 120061 ) . except for the value of the coefficient. 



The dynamical stability of the surface density profile given by equa tion (jll) must be 
checked using the well-known Rayleigh criterion (e.g.. IChandrasekhaiill96ll ). For the surface 
density profile to be stable, in the local coordinate system. 



dx - p 



(7) 



must be fulfilled. The velocity Vy can be obtained by solving the equation of motion for the 
X component. 



3Qlx + 2npVy - c 



dx ' 
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where c is sound speed. Substituting equation (jl]) into equation ([8]), we obtain 



X 



where h is the disk scale height defined as h = c/VL^. From equations ([7]) and 
that condition ([7]) is not fulfilled for x < Xm, where 



121/5 



0.207 



2/5 



121/5 



h 



3/5 



O.lr, 



2/5 



10-V21]p 



-1/5 



10-3M, 



2/5 



P- 



(9) 
we find 

(10) 
(11) 



In the region x < x^, the density gradient thus has to be small enough to satisfy the 
Rayleigh condition: the density profile would be relaxed to be marginally stable for the 
Rayleigh condition. Integrating dvy/dx = —2Qp inward from Xm we obtain 



Vy ^y,m 217p(x ^m) 

in the region x < x^a- Substituting this into equation ([8]), we finally obtain 

^(") = ^-^"p1^-2U~4tJ +32 It) 

for X < Xm- 

In summary, the equilibrium profile of surface density is given as 

Sr(x) for X < Xm, 



(12) 



(13) 



^(^) \ ( \ ? ^ 

y Lvis(a;) for x > Xm- 

An example of the density profile is shown in Fig{TJ 



(14) 



As mentioned above, we have assumed Vy = — (3/2)f2pX in deriving equation (jlj). How- 
ever, based on equation ([9]), the velocity at 



(15) 



This means that the maximum deviation from the Keplerian shear velocity, — (3/2)f2pX, is 
less than 10% of the absolute value. Hence the assumption of the Keplerian shear velocity 
is reasonable for x > Xm- 



- 8 - 



Finally we briefly describe the dependence of the surface density on the parameters. 
Figure [2] shows the surface density as a function of x/r^ for several different values of plan- 
etary mass Mp (a), viscosity v (b), and scale height h (c). Figures [2^ and 2b show the gap 
becomes deeper and wider with increasing planetary mass and decreasing viscosity. The 
width and depth of the gap are determined by competition between gravitational scattering 
by the planet and viscous diffusion of the disk gas. The gap becomes wider and deeper when 
planet's gravity is stronger or viscous diffusion is less efficient. From a mathematical point 
of view, this is because both S^is and depends on Mp and u in the form of (Mp/z/)^/'^ 
through I. As for the dependence of S on /i (see Fig. [21:;), all the curves are independent of 
h for X > Xjn, where S = Svis which is independent of h (see eq.[l]). In other words, surface 
density in a; > Xm is determined by the balance between viscous torque and gravitational 
torque from the planet and has nothing to do with pressure. On the other hand, difference 
in surface density appears at x < x^, where S is determined by the Rayleigh condition (see 
eq.[T3]) and thus depends on all of h, z/, and Mp through Xm- 



3. Gas accretion rate onto planets 

In this study, we describe the gas accretion rate onto the planet, Mp. Through a 
series of local high-resolution hydrodynamic simulations of the gas accretion flow onto giant 
planets in protoplanetary disks, TW02 obtained an empirical relation between accretion rate 
normalized by unperturbed surface density and /i/rn, which is the only one parameter of the 
local system (see eq. [18] of TW02). Furthermore, TW02 demonstrated that only the gas in 
the two bands at |x| ~ 2rH accretes onto the planet. Thus, writing equation (18) of TW02 
in the explicit form, for 0.5 < h/r^ < 1.8, we have 

Mp = is,,,, (16) 

where 



h\~' /M 



4/3 

which is the area, the gas in which is to be accreted onto the planet per unit time (hereafter 
the accretion area), and Sacc = S(2rH). 

The reason why small h or large Mp yields high Mp is explained as follows. Disk gas 
loses its energy by passing through spiral shocks around the planet and consequently accretes 
onto the planet. Small h (i.e., small sound velocity) or large Mp (i.e., strong perturbation 
on gas motion) yields a large Mach number of the disk gas, which results in strong shocks. 
In addition, an increase in Mp expands the planetary feeding zone, resulting in large Mp. 



- 9 - 



Note that the position of the accretion band, Xacc, is located at ~ 2rH in the highest- 
mass case of TW02 (i.e., Ciso = 0.5 in their notation; see Fig. 7 of TW02 where Xacc is 
shown to depend shghtly on h/rn). When Ciso < 0.5, a gap exists around the planetary 
orbit. In such situations, appropriate choice of Xg^cc is crucial in determining Mp because a 
small difference in x^cc makes a large difference in Eacc and thus Mp, while choice of Xacc has 
little influence on Mp in low-mass cases. That is why we adopt the value of Xacc(= Srn) for 
the highest-mass case of TW02. 



Note also that, as shown by iD'Angelo. Kley. fc HenningI (l2003l ). the gas accretion rate 



obtained by two-dimensional simulations are usually higher than that obtained by three- 
dimensional simulations. Thus the accretion rate given by equation (|T6|) could be overesti- 
mated. We will discuss this issue in §4.4[ 



4. Evolution 

In this section we show the evolution of the growth rate and mass of an accreting planet. 
To gain a proper understanding of the basic behavior of the planetary accretion, we first 
explore two simple cases with no disk dissipation (i.e. constant Soo) and with exponentially- 
decreasing in § §4.11 and 14.21 respectively. Then, we investigate the accretion of a giant 
planet embedded in a viscously evolving protoplanetary disk in §4.31 

The numerical procedure is as follows. Except in §4.H we first calculate the unperturbed 
surface density, Eqo, from equation (12T!) in §4.2l or from equation (126!) in §4.31 For given values 
of the parameters, h/r^ and u/^r'^Qp), we then calculate Eacc and A, using equation (jlj) if 
a^acc > OY cquatiou (IT^ if Xacc < a^m- Finally we integrate Mp with respect to time using 
Sacc and A (see eq. [16]) to obtain the time evolution of the planetary mass. 



4.1. Case without disk-dissipation 

We first show the evolution of the planetary mass without global disk dissipation. In 
this case the input parameters are SooTp/M^,, h/vp, and u/^r'^flp). 



4.1.1. General Properties 



Figure [3] shows the evolution of the planetary mass and the gas accretion rate for several 
values of the three parameters. Without global disk depletion, the only way to suppress gas 
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accretion is opening a deep gap around the planetary orbit by the planet itself. As seen in 
Fig. [3], the evolution can be divided into two phases. The first phase is growth without a 
gap (pre-gap phase), while the second phase is growth with a deep gap (post-gap phase). 
Properties of the evolution in both phases can be explained in analytical ways below. 

In the pre-gap phase, Mp is almost proportional to Mp^ (see Figs. [3]-(c),(f),(i)). In this 
phase, the planetary mass is insufficient to open a gap, so that Sacc — '^oo- Then we can 
easily integrate equation (ITBi) with equation (ITTI) : 



h 



p,init 



-1/3 



S 

3" 



-2 



4 ' 



where S corresponds to 0.29 in equation (ITTI) and Mp init is the initial mass of the planet that 
corresponds to the mass beyond which the gas accretion is limited by disk-gas supply. This 
equation implies that if Sacc is constant, the gas accretion rate diverges at a time 



3 
5 



M, 



p.init 



4.8 X 10^ 



-1/3 



M 



p,init 



1/3 



-1 



^ 1 



10-5M, 



10 



10 



-1.5^ 



-1. 



(19) 



Tdiv corresponds to the end of the pre-gap phase. 

In the post-gap phase, Mp decreases with Mp in an exponential fashion (see Figs. [3]- 
(c),(f),(i)). Such a steep decrease in Mp is due to a steep decrease in Sacc with respect to 
Mp via £, Xm, and rn (see eqs.|l] and [IS])- We can derive approximate expressions for the 
accretion rate in this phase, which shows Mp is inversely proportional to time as shown in 
Figs. [3](a)-(c) (see appendix for the detail). 



4.1.2. Dependence on the parameters 

We now see the dependence of the evolution of planetary mass on the unperturbed 
surface density, Sqo (Figs. [3]-(a),(b),(c)). The evolution timescale is inversely proportional 
to Soo (e.g.. Fig. [3]-(b)). This is simply because Sacc is proportional to Eoo, whereas A is 
independent of Soo- The peak value of Mp is thus proportional to Eoo accordingly (e.g.. Fig. 
[1(c)). 
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The dependence on scale height, h, is shown in Figs. [3]-(d),(e),(f). In the pre-gap phase, 
Mp is simply proportional to /i"^ and the evolution timescale is thus oc h"^. This is because 
the accretion area A oc (see eq.[T7]) in the pre-gap phase where Sacc = ^oo- Even in the 
post-gap phase, the dependence of Mp on h is the same as that in the pre-gap phase if h 
(i.e., Xm) is small {h/r < 0.03). As long as h is so small that Xm < 2rH, Eacc is determined 
by Evis and is thus independent of h (see eqs.jl] and [H]). When h is so large that Xm > 2rH, 
Sacc = Sr, which depends on h. Hence Sacc increases with h for large h. The transition 
between the low-h and high-/i cases takes place when Xm = 2rH, namely, from equation (II ip : 

The dependence on viscosity, z/, is shown in Figj3l-(g)(h)(i). In the pre-gap phase, the 
evolution of Mp is almost the same for different values of because Eacc is almost constant 
(^ Soo) in the pre-gap phase and A is independent of v. Thus the timing when Mp reaches 
at the peak value does not depend on v either (see eq.[l9]). In the post-gap phase, Mp at a 
given t increases with u. In appendix [HI we have derived approximate solutions for Mp in the 
post-gap phase. From equations flB2p and flB3p . one finds that Mp oc in high-z/ cases (i.e., 
a^acc > a^m)- In low-z/ cases (i.e., Xacc < a;m), Mp also increases with i/, but the dependence 
is rather weak (see eqs. |B7] and |B8j ). Also, the accretion rate is inversely proportional to 
time (see eqs. |B3] or [B8] ). so that the maximum accretion rate can be roughly estimated 
by eqs. |B3j or |B8j aX t = Tdiv 



4.2. Case with exponential disk dissipation 

Next we examine the evolution of the planetary mass in a simple case where the disk 
surface density decreases in an exponential fashion with a time constant of r^cp'- 

^oo = Soo.init 6Xp I 1 , (21) 

where init is the initial surface density at infinity. 

Figure m shows the evolution of Mp and Mp for Xdep = lO^fi" ^ and three values of T^oo^mit- 
In the high-Soo,init case (Eoo,init/(^*'^p ^) = 10""^-^), Mp increases with time to reach a peak 
at t ~ 1.5 X 10^f2p^ and then decreases with time, which is similar to the evolution without 
disk dissipation shown in §4.11 This is because Tdiv < Tdep in this case. In the low-Soo,init 
case (Soo,init/(A^*'"p ^) = 10~^'^), on the other hand, Mp decreases without experiencing a 
significant increase that occurs in the high-Soo,init case. This is because Tdiv > Tdcp in this case 
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(rdiv oc ^^init'^ see eq. [19]). Because of such different evolution of Mp, the mass evolution also 
differs between high- and low-Soo,init cases (see Fig. HJd): The planet captures a significant 
amount of gas to be a Jupiter-like planet in the high-Ego, init case, whereas the planet captures 
only a small amount of gas to be a Neptune-like planet in the low-Soo,init case. The boundary 
between the two regimes is determined by the condition T^ep = Tdiv, namely, from equation 

Figure \5\ shows the final mass defined by 

poo 

Mfinai = I M^dt (23) 
as a function of a quantity ( defined by 







^oo,init \ / ^dep 



^ - [km) w) ■ '''' 

The dependence of Mgnai on r^cp is the same as that on Soo,init because the ratio of rjiv 
(oc Sj^^jjjjJ to Tdep determines the evolution. Figure [5] illustrates that when ( is small (< 1), 
the final mass is almost the same as the initial mass. In this regime, the planet captures 
only a tiny amount of gas and becomes a Neptune-like planet. When ( > 10, on the other 
hand, the planet captures substantial gas to be a gas giant planet like Jupiter and the final 
mass does not depend on the initial planetary mass. This is because the planet becomes 
large enough to form a gap and the evolution slows down significantly after gap formation; 
the condition for gap opening is determined by the planetary mass. The transition occurs 
when Tdiv = Tdep, namely, from equation f|T9|) . 

C..4.8feV"7^V. (25) 



Figure [6] shows that the final mass as a function of ( and u in the cases with h/r = 0.1 
and h/r = 0.032. A portion characterized by a steep gradient is illustrated to shift to the 
left with decreasing h/rp {(t ~ 5 for h/vp = 0.1, while Ct ~ 0.1 for h/vp = 0.032, see eq. |25j). 
On the left side to the portion (low-^ case) , the final mass is almost equal to the initial mass 
and thus almost independent of viscosity. On the right side, the final mass increases with 
viscosity. This is because the growth is limited by opening a deep gap; the mass for gap 
opening depends on viscosity. 
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4.3. Case with global disk evolution 

4-3.1. the model 

We have so far neglected global evolution of the disk, although we examined a simple 
case with an exponential depletion of the disk gas in §4.21 In reality, the unperturbed surface 
density, Sqo, changes because of global disk evolution. In this section, we examine the effects 
of global disk evolution on planetary growth. 

Because a detailed description of disk evolution is beyond the scope of this paper, we 
assume that the disk surface density changes with time in such a way 

Eoo = Sss(r,t)exp I ^ j , (26) 

instead of equation (I2T1) : represents a change in the surface density due to viscous 
diffusion of the disk gas and e~*/'^'^'=p is introduced to mimic photoevaporation of the disk. 
For Ess, we adopt the self-similar solution with a prescription for disk viscosity given by 



(IHartmann et al.lll998h 



Sss(r,t) = -^ C/^exp , (27) 



^disk / ^\ 3/2 (_r_/Ro 

where M^isk is the initial disk mass and i?out is the initial disk size; fgs is defined as 

Tss = — + 1, (2^ 

''Vis 

with a typical timescale of global viscous evolution 

p2 

_ ^ont 



3t'out 



^5 ( " A ^ i ^lAU \ / -Rout 



where z/out is viscous coefficient of the disk gas at r = -Rout- In deriving equation (127|1 . we 
have assumed a temperature distribution, T oc r"^/^, which results in oc r. 

Using the equations above, we calculate the local gas accretion rate onto the planet, 
denoted hereafter by Mp^iocai, in the same way as we did in §4.2[ However, when Mp^iocai 
is larger than the radial mass flux due to viscous diffusion in the disk, the latter limits the 
planetary growth. The mass flux through a ring with radius r is given by 



Mdisk(r-,t) = 27rrSss(r,t)fr 



2 ^ss-Rout / V ^ss-Rout 



^sl'/'expl-— — ), (30) 
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with the radial drift velocity of the diffusing gas in a Keplerian disk 

3z/ /SlnfSz/) r 



Vr 



r \ dlnr 

Thus we calculate the gas accretion rate onto the planet in such a way 



^p,local if ^p.local < ^disk 



disk 



if M, 



p, local 



> M, 



disk' 



(31) 



(32) 



There are several studies that modeled gas accretion rate onto gas giant planets. iGuillot fc Hueso 



( I2OO6I ) used a formula in which the gas accretion rate onto the planet is 0.3 times the radial 
mass flux of the disk gas due to global visco us diffusion. However their formula is an em- 
pirical one and apphes in a limited situation. iLubow. Seibert. fc Artymowica (119991 ) carried 
out a series of hydrodynamic simulations to obtain gas accretion rate onto planets. Their 
simulations demonstrated that the planetary accretion rate can be larger than the diffusion 
flux because of gradients imposed by the gap. However, what they found is probably a 
transient phenome i ion th at occurs on a timescale much shorter than the viscous timescale. 
Veras fc Armitaed (120041) used an appro x imate forn aula based on hydrodynamic s imula tions 
done by lLubow. Seibert. fc Artymowica (119991 ) and lD'Angelo. Henning. &: Kleyl (120021 ). 



4-3.2. dependence on the disk parameters 

Figure [7] is a similar figure to Fig. [6] and shows the final mass of the planet, Mgnai, as a 
function of the initial surface density, Eoo,init, and viscosity, z/, in the case where h/r = 0.1 
(left panel) and 0.032 (right panel). Rout = lOrp, and r^ep = lO^fip^. The value of r^ep 
corresponds to C = 10^(Soo,init/(A^*'"p ^))5 so that the ranges of this figure are the same as 
those of Figure M for the horizontal axis as well as the vertical axis. 

Since we have adopted the self-similar solution for the surface density evolution of the 
global disk, there appear two additional limits to the final mass. One arises from the ap- 
pearance of the viscous timescale r^is, in addition to the timescale of exponential decay rdcp. 
When Tvis < Tdep, the effective disk lifetime is r^is- Thus, in order for planets to be massive, 
the growth timescale Xdiv should be shorter than Tvis, so that the additional limit is given by 

''Vis ''"div 

r2(]p \l0-mj \0.1rJ \lO-m,/rl) \l0rj ' ^ ^ 

which runs from the lower left to the upper right in Fig|7| because u is linearly proportional 
to the surface density. As described above, the effective disk lifetime, rufetime, is the shorter 
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one of the two timescales, so that we here define 




(34) 



The other hmit is the total mass of the disk. This hmit does not appear in FigjG] (i.e., the 
case without global viscous evolution). This can be clearly seen in the case of h/r = 0.032 
on the right panel of Fig. [7] where contour lines are vertical. Note that the steep boundary 
(being equivalent to C = Ct) characterized by r^iv = Tdep that is seen in Figj6]is not found 
in Fig. [3 This is simply because the two limits described above are more severe than the 
condition of rjiv = Tdep in the cases shown in Figure [3 However, when Tdep is smaller, the 
boundary determined by r^iv = Tdep moves rightward, so that the vertical steep boundary 
determined by Tdiv = Tdep emerges even in the case with global disk evolution described in 
this subsection. 



In Figure Ui the dependence of the final mass was shown for the normalized disk param- 
eters. However, one might want to know its dependence on semi-major axis. We thus show 
the final mass as a function of semi-major axis in Figure [S] for a typical case with a = 0.01, 
h/r = 0.032 at lAU, Tdep = lO^r, ^out = lOOAU, Mp^i^it = 3.2 x lO'^M,, which corresponds 
to 10 Earth masses in the solar mass system, and Mdisk = 0.013M*, which yields a surface 
density roughly 0.1 times that of minimum mass =olar nebula at lAlfl. In addition to the 
final mass given by equation (!23l) . two kinds of the mass defined by 



are shown to clarify the limiting process for the final mass. 

As shown in Figure [8], the final mass (solid line) is determined basically by the smaller 
of Mpjocai (dotted line) and Mp^disk (dashed line), although the limiting process varies with 
time. Thus one finds three characteristic regions that are divided at the intersection of the 
two curves. 

In the inner region (rp < lAU), Mfinai is almost equal to Mp^iocai and is proportional to 

1 /2 

Tp . The duration of the pre-gap phase is shorter than the disk lifetime (i.e., Tdiv < Tufetime) 



4-3.3. classification by semi-major axis 




and 




(35) 



^This surface density is much lower than the minimum mass solar nebula model because the slope of the 
surface density is shallower and the disk size is larger. 
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in this region. The gas accretion is thus suppressed because of gap opening. This region 
is hereafter called the gap-limiting region. In this situation the final mass can be roughly 
estimated by integrating Mp from Tdiv to riifetime- Since Tdiv < Tufctime in the gap-limiting 
region for the values of the parameters used in drawing Fig. [HI Mfinai is readily calculated as 
(see Appendix [B]) 

-'Tdiv 

= Mtrans log ( Tl^^^^ j = Mfinal.gap (36) 

where Mtrans; defined by equation (lB2p . is a planetary mass when Sacc = ^00/ e and 2rH < x^a- 
Except for the weak dependence on the other parameters included in the log term, the final 
mass given by equation fl36|) is determined principally by u and increases with u. This is 
because larger planetary mass is required for the planet to open a gap in the case of higher 
viscosity. Since Mtrans oc ///(r^fip) and oc in our disk model, Mtrans oc arp'^^, which is 
consistent with the result shown in Fig. [HI 

In the intermediate region (1AU< Vp < lOOAU), Figure [8] illustrates 

Mfinal ^ Mp_disk (37) 

In the region, viscous diffusion of the disk gas limits the gas accretion onto the planet. This 
region is hereafter called the diffusion-limiting region. In this region, radial mass flux due to 
viscous diffusion is uniform where Vp <C -Rout, so that the final mass is usually insensitive to 
Tp. In the case with r^is < r^cp, which is the case of Fig. [HI most of the disk gas is transferred 
toward the central star by viscous diffusion, thus the final mass is roughly the disk mass 
(~ 0.013M*). Hence we have 

Mfinal.diff ~ Mdisk- (38) 

On the other hand, when r^i^ ^ Tdep, evaporation dominates viscous diffusion in disk dissi- 
pation. In this situation, the final mass at r <^ i?out is given by 

/"Tdcp 

Mfinal,difT ^ / I Mdisk 

^ , , Mdisk. (39) 

Note that, even when Tvis ^ r^ep, the gas accretion rate onto the planet is limited by viscous 
diffusion. 

In the outer region (r^ > lOOAU), the final mass suddenly decreases to Mp^init. This is 
because the planetary growth finishes in the pre-gap phase because of the dissipation of the 
disk gas (i.e., Tdiv ^ Tufetime)- We call this region the no-growth region. 







1 / Tdep 
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The positions of the boundaries between the three regions are given as follows. At 
the boundary between the gap-limiting region and the diffusion-limiting region (its position 
being denoted by rb), Mgnai.gap = ^finai,difr- When r^is < Tdep, from equations and ([3SD, 
we obtain 

~ 5 ' ( _J}m \ ' ( ^disk \ ' / log(rvis/rdiv) \ 

"-Hl0-2y' *^10-i-5AUy' \ h ) ' ^ ' 

Here we have assumed the log factor in equation fl36|) is constant. When TvIs > Tdep, from 
equations ( l36l) and ( l39i) . we have 



^ 1 « A~V ^lAU ^ V Mdisk ^ VrdepV / log(rdep/rdiv) ^^ ^ 

The location of the boundary between the diffusion- limiting region and the no-growth region 
is basically determined by Xdiv ~ Tufctimc and is denoted by r^. When r^is < Tdep, from 
Tdiv = Tvis, we obtain 

1/3 / V ^ ... \ / h ... \ -4 



e — 



70 ( ( ^"o.lAU \ ( /^lAU y ( Ro-t \ (AO) 



When Tvis > Tdep, from Tdiv = Tdep, we obtain 



1/3 / V. ^ ... X / , ... X -2 



^^°(,io-5mJ (io-w^/TauO (lo-'^^Au) (lov)"^^- ^^^^ 

Figure [9] shows the final mass for several different values of the disk parameters. The 
discussion above suffices to understand this figure. Note that there is a hollow at ~ 50 AU 
in the case where the value of a is lower by a factor of 10 than its nominal value (thick 
short-dashed line). This is, however, artificial. In the self-similar solution we have adopted, 
the direction of the radial gas flow due to viscous diffusion changes at r = (_Rout/2)fss, which 
means that the mass flux is zero at the point, and the accretion rate onto the planet is also 
zero accordingly (see eq. [32]). Since f^s increases with time (see eq.[2B]), the point moves 
outward with time. When Tdep > t^is, the point moves significantly with time and integration 
of mass flux at a particular point eliminates the effect. In contrast, when Tdep ^ tvis, the 
point hardly changes before the exponential depletion; the planet around r = -Rout/2 thus 
cannot grow significantly. 

Another important aspect is on the final mass in the gap-limiting region. We can see 
that the rp depe ndence of Mfjnai agrees with th at of the mass given by the viscous condition 



for gap opening (jLin &: Papaloizoulll979l . Il993l ): 



AW = «(^)Af. (44) 
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The dependence also agrees with that of Mtrans, which is the mass when Sacc is (l/e)Soo in 
the case of 2rH > (see eq. |B2j ). and re- written as 

- 9 X 10- (^) [-^Y" M.. (46) 

Both of the two masses indicate gap-forming masses and they agree with each other within a 
factor of about 2. However the final mass shown in our model is larger by a factor of 10 than 
^p,vis or Mtrans- This is because the two masses correspond to the masses at which a gap is 
about to form and does not mean the masses at which growth stops. A further increase in 
Mp after reaching Mp^vis or M^ans is reflected in the log term in equation ( l36l) . 



4.4. Remarks on the uncertainties in the model 



We have assumed that the distribution of the disk gas around the planetary orbit (i.e., 
the surface density) is always in the equilibrium state that is determined by the balance be- 
tween viscous stress and gravitational scattering by the planet. However, for the equilibrium 
state to be achieved, non-equihbrium distributions must be relaxed by diffusion. The diffu- 
sion timescale, Tvis.iocai, is estimated by the gap width (~ 2 x 2rH) divided by the viscosity 
as 



T-vis,Iocal = 7.7 X 10' 



10-3M, 



2/3 



lo-v^^p 



-1 



-1 



(47) 



This timescale should be compared with the typical growth timescale (i.e., Tdiv; see eq. [E]). 
The comparison indicates that Tvisjocai < ^div in niost of the cases presented in this paper. 
Hence the assumption is appropriate unless quite low viscosities (yielding large Tvisjocai) or 
high surface densities (yielding small taw) are considered. 

The empirical formula for the accretion rate based on local two-dimensional isothermal 
hydrodynamic simulations (TW02) can be different from that derive d based on more realistic 

simulations including, for example, three-dimensional accretion flow (iD'Angelo. Kiev. &: Henning 

2OO3I : iBate et al.ll2003l ). non- isothermal equation of state (TW02), and magnetic field (IMachida. Inutsuka. & 
2006h . The modification would change the form of rdiv, which would yield quantitatively dif- 
ferent results. However, even in that case, following our prescription given in this paper, 
one can easily calculate the mass evolution of planets in a similar way. The formula given 
by TW02 can be considered as the highest limit; thus, for example, is expected to shift 
inward if more realistic models for the accretion rate are used. We compare the accretio n 
rate eq. [16] with those obtained by global simulations (lKleylll999l : iLubow &: D'Angeldl2006l ). 
Although the value of our accretion rate are usually larger (by up to a factor of 10) than 
those given by the simulations, the dependence on viscosity and planetary mass is consistent. 
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We have not taken into account a limit in terms of geometry of the accretion flow, 
namely a within which equation f|T71) is applicable. We can make a simple estimate of the 
maximum accretion rate if we assume all the gas within < x < 2^/3ru in y > and 
<x<Omy<0 accrete to the planet, where x = 2-\/3rH at ?/ rn is the point 
where potential energy on the rotating frame is the same as that at the Lagrange points LI 
or L2. In this case, the maximum accretion rate (normalized by surface density) is given by 
3^/^6rpfip(Mp/M*)^/^ and thus A (eq.[l7]) becomes larger than the maximum value when 



Mp > 5.6 X lO-^M, =Mp,ent. (48) 



Although planetary mass can be larger than Mp ^rit depending on the parameters, such a 
massive planet should already has a deep gap, so that the accretion rate is reduced greatly 
when the planetary mass reaches Mp^crit- Consequently the geometric effect has no signif- 
icant influence on the final mass. Note that the discussion here is based on an approx- 
imation that streamlines can be well described by particle motions on the framework of 
the restricted three-body probleni . Since this approximation is not valid when r-^ < h 



( iMasset. D'Angelo. and Kleyll2006[ ). it is not applicable for small mass planet's cases. How- 
ever, we consider the upper limit of applicable mass, where the planet mass is, in most cases, 
large enough to satisfy rn > h, so that the approximation is valid for the situations that we 
consider here. 

In this paper, we have put off the issue of planetary migration. Since we consider a phase 
in which a gap exists around the p lanet, we may have to consider the type II migration (e.g. 



id) 



WardI 119971 : iWard and HarnI 120001 ). especially in the gap-limiting region (e.g., < 1 AU in 



Figure [8]). The deep gap created by the planet blocks the accretion flow toward the central 
star. It follows that the planet is pushed inward by the gas exterior to the planet's orbit. 
Inclusion of the effect of planetary migration is our future work. 



5. Summary 

To gain a systematic understanding of the final masses of gas giant planets, we have 
simulated the long-term accretion of gas giant planets after the onset of the supercritical gas 
accretion in a variety of situations, depending on four disk parameters such as disk mass, 
viscosity, scale height, and semi-major axis. To do so, we have made a semi-analytical model 
to simulate the mass evolution, which enabled us to study the final mass of gas giant planets 
for extensive ranges of all the parameters. 

We have first made a ID analytical model of the equilibrium surface density profile 
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around a protoplanet, from consideration of the balance of torque and the dynamical stabil- 
ity O- Combining the surface density profile with an empirical formula for gas accretion 
rate that was obtained on the basis of hydrodynamic simulations by TW02, we have obtained 
a formula for gas accretion rate as a function of planetary mass, viscosity, scale height, and 
unperturbed surface density (^. We have then integrated the gas accretion rate numeri- 
cally with respect to time to simulate the long-term accretion of gas giant planets (§4). To 
understand the basic behavior of the planetary accretion, we have explored two simple cases 
with no disk dissipation ( §4.ip and with exponentially-decreasing surface density ( §4.2p . Fi- 
nally, we have simulated the long-term accretion of gas giant planets embedded in a viscously 
evolving and evaporating disk to obtain the final mass of gas giant planets as a function of 
semi-major axis ( §4.3p . We have consequently found the following three different regions 
depending on limiting processes on the final mass. 

In the inner region (rp < ry,; see eqs. |l0] & [S]), the planet grows quickly to form a 
deep gap to suppress the gas accretion from the disk by itself within the disk lifetime {gap- 
limiting region). We have found the final mass in this regio n is roughly 10 times la rger than 



that determined by the viscous condition for gap opening (iLin &: Papaloizoulll993l ). This is 
because the condition for gap opening only expresses the condition when a gap begins to 
form, and is by no means equivalent to the condition that the growth is terminated. 

In the intermediate region (rb ^ rp < see eqs. |12] & [43]), radial transfer of the disk 
gas toward the planetary orbit limits the gas accretion before the planet opens a deep gap; 
the final mass is thus limited by viscous diffusion of the disk [diffusion-limiting region). We 
have found that when the evaporation timescale r^ep is shorter than the viscous- diffusion 
timescale r^is, the relationship between the final mass Mgnai and the disk mass Mdisk is given 
by Mfinai ~ (l/2)(rdop/rvis)Mdisk, whereas Mfi^ai ~ Mjisk when r^ep > r^is. 

In the outer region (r > re), the planet captures only a tiny amount of gas by the time 
the disk gas completely dissipates {no-growth region). Saturn and possibly Uranus/Neptune 
are likely to have experienced the situation. 

In this study, we have gained a clear understanding of the final masses of gas giant 
planets, deriving analytical expressions for them in three characteristic regions (eqs. [5B]-|39j) 
and the locations of the boundaries between the three regions (eqs. |10]-[l3]). To understand 
the mass-period distribution of gas giant planets in extrasolar systems found by radial- 
velocimetry, we need to take several additional processes into consideration. Especially, 
planets in the gap-limiting region would be susceptible to the type II migration, because the 
gas exterior to the planetary orbit blocked by the gap pushes the planet inward. Inclusion 
of planetary migration is our future work. Also, inclusion of core accretion processes and 
the gas accretion process governed by the Kelvin-Helmholtz contraction of the envelope is 
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needed especially to determine the initial mass and the origin of time of our model. 



Although the final masses of gas giant planets were focused on in this paper, the ac- 
cretion process for reaching the final mass is also important to resolve issues relevant to 
planet formation. Growing giant planets dynamically affect other bodies in a planetary 
system. Satellites are likel y to form in sub-disks around accreting gas giant planets (e.g. 



Canup fc Wardll2002l . I2OO6I ). The long-term accretion of gas giant planets may affect the in- 



ternal structure and evolution of isolated young gas giants (IMarley et al.l 120071 ). This would 
be important for future direct detection of young gas giants. 
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(MEXT), Grand-in-Aid for Scientific Research on Priority Areas, "Development of Extra- 
solar Planetary Science" (MEXT-16077202). 



A. Disk properties 

Typical values of the disk parameters in our model are summarized with their depen- 
dence on the semi-major axis of the planet. We assume that the disk is optically thin and 
i n the radiatiy e equilibrium with the stellar radiation. The disk temperature is given by 



flHavashil 119811 1 



where and Lq is luminosities of the central star and the sun, respectively. Hence the 
normalized scale height (i.e., the aspect ratio of the disk) is 

- = 0.033 (A2) 



rp VlAU 

Note that the ratio is independent of the stellar mass because we assume the stellar luminosity 
is proportional to the fourth power of stellar mass, which is observationally known in the 
case of main sequence stars. 



On the a-prescription {a = u/jch)), which is widel y used to describe poorly- known tur- 
bulent viscosity in the disks (jShakura fc Sunyaevlll973l ). the normalized viscosity coefficient 
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with the temperature distribution given by equation flAll) is 



If we assume that disk surface density is proportional to the mass of the central star, nor- 
malized surface density is 

1.9x 10-^f^)'^Y^V (A4) 



M,/r2 V1AU7 Vl 

where / is a factor relative to the surface density of the minimum mass solar nebula model 
when the mass of the central star equals to that of the Sun. The dissipation timescale of the 
disk is 

3/2 f M,\ "^/^ , , , 

(A5) 



^dep 9^ V in6 ^dcp ^ / r-p \ -3/2 / M, \ 



where Mq is the stellar mass. The accretion rate is 



4.8 X 10"^ L^_o J / T7^ (t^) (A6) 



B. Approximate solutions for the accretion rate in the post-gap phase 

The surface density of the gas in the accretion band at x ~ 2rH, Sacc, is given by S^is, 
when the accretion band exists in the region with the dynamically stable profile of surface 
density determined by the balance between viscous diffusion and scattering due to the planet, 
namely, when 2rH > x^- Substituting x = 2rH into equation (jlj), one obtains 

Sacc = Svis(2rH) = Soocxp ( - J , (Bl) 

V trans / 

where 

Mtrans = 2771 (^-^^ M,. (B2) 

Using equation (IBip . one can rewrite Mp = Sacc^ as 

Mp ~ (B3) 
for large t and Mp, if one neglects the weak dependence of A on Mp. 
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When 2rH < Xm, Mp depends on Mp in a more complicated manner. As seen in equation 
f|T3|) . S/j is not a simple exponential function of Mp unlike S^is. Defining a function f{Mp) 
as 

one can write Sacc in the form 

S,,e = SR(2rH) = Sooexp(/(Mp)). (B5) 

Since Mp does not change so much in the post-gap phase, one only has to integrate Mp in 
a limited range of Mp. Referring to AMp as an increment from a mass, Mp^ (i.e., Mp = 
Mp o + AMp), one expands e^^ in terms of the small quantity AMp in such a way 

Eacc = Sr(Mp,o + AMp) = Sooexp (^/(Mp,o) + ^AMp^ 



where 



c/Mp Mp 



Thus one obtains 

i-df/dMp)-^ 

Mp ~ ^ , (B8) 

V 

for large t and Mp. 



REFERENCES 

Alibert, Y., Mordasini, C, & Benz, W. 2004, A&A, 417, L25 
Artymowicz, P., & Lubow, S. H. 1996, ApJ, 467, 77 

Bate, M. R., Lubow, S. H., Ogilvie, G. I., & Miller, K. A. 2003, MNRAS, 341, 213 

Bodenheimer, P., & Pollack, J. B. 1986, Icarus, 67, 391 

Boss, A. P. 2002, Earth and Planetary Science Letters, 202, 513 

Boss, A. P. 1989, ApJ, 345, 554 

Cameron, A. G. W. 1978, Moon and the Planets, 18, 5 



-24- 



Canup, R. M., & Ward, W. R. 2002, AJ, 124, 3404 
Canup, R. M., & Ward, W. R. 2006, Nature, 441, 834 

Chandrasekhar, S. 1961, Hydrodynamic and Hydromagnetic Stability, Clarendon press 
Crida, A., Morbidelli, A., & Masset, F. 2006, Icarus, 181, 587 
D'Angelo, G., Henning, T., & Kley W. 2002, A&A, 385, 647 
D'Angelo, G., Kley, W., & Henning, T. 2003, ApJ, 586, 540 
Guillot, T., & Hueso, R. 2006, MNRAS, 367, 47 

Hartmann, L., Calvet, N., GuUbring, E., & D'Alessio, R 1998, ApJ, 495, 385 

Hayashi, C. 1981, Prog. Theor. Phys. SuppL, 70, 35 

Ida, S., & Lin, D. N. C. 2004, ApJ, 604, 388 

Ikoma, M., & Genda, H. 2006, ApJ, 650, 1150 

Ikoma, M., Nakazawa, K., & Emori, H. 2000, ApJ, 537, 1013 

Kley, W. 1999, MNRAS, 303, 696 

Lecar, M. & Sasselov, D. D. 2003, ApJ, 99 

Lin, D. N. C., & Papaloizou, J. C. B. 1979, MNRAS, 186, 799 

Lin, D. N. C., & Papaloizou, J. C. B. 1993, Protostars and planets III, 749 

Lubow, S., H., Seibert, M., & Artymowicz, P. 1999, ApJ, 526, 1001 

Lubow, S., H. & D'Angelo, G. 2006, ApJ, 641, 526 

Machida, M. N., Inutsuka, S., & Matsumoto, T. 2006, ApJ, 649, 129 

Marley, M. S., Fortney, J. J., Hubickyj, O., Bodenheimer, P., & Lissauer, J. J. 2007, ApJ, 
655, 541 

Masset, F., D'Angelo, G., & Kley, W. 2006, ApJ, 652, 730 
Mizuno, H. 1980, Prog. Theor. Phys., 64, 544 

Pollack, J. B., Hubickyj, O., Bodenheimer, P., Lissauer, J. J., Podolak, M., & Greenzweig, 
Y. 1996, Icarus, 124, 62 



-25- 

Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 
Tanigawa, T., & Watanabe, S. 2002, ApJ, 586, 506 (TW02) 
Veras, D., & Armitage, P. J. 2004, MNRAS, 347, 613 
Ward, W. R. 1997, Icarus, 126, 261 

Ward, W. R. & Harn, J. M. 2000, Protostars and planets IV, 1135 



This preprint was prepared with the A AS IM^jX macros v5.2. 



- 26 - 




1 2 3 4 5 
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Fig. 1.— Surface density profile for Mp/M* = lO^^^ v/{rl^v) = 10"^ and h/r^ = 0.1. Solid 
line is S^is given by equation (jl]) and dashed line is given by equation f|T3l) . The two 
curves connect smoothly with each other at 2.08h in this case). Thick line 

is the actual surface density we use in this paper. 
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Fig. 2. — Surface density profiles as a function of distance from the planetary or- 
bit. Top panel (a) shows planetary mass dependence (Mp/M* = 10~^ (solid line), 
10^'^"'', 10^'^, 10^^'^, 10^^ (dot-dashed line)), middle panel (b) shows viscosity parameter de- 
pendence (z//(rpr2p) = 10^^(solid line), 10^^, 10^''', 10^^(dotted line)), and bottom panel (c) 
shows scale height dependence {h/r-^ = 10^^(solid line), 10~^'^, 10^^, 10^°-^ (dotted line)). 
The standard values of the parameters are Mp/M^ ~ 10~^, v/{rpVL.p) — 10~^, and 
h/r-p — 10~^. The position of the symbols -|- indicates where the top panel and 

the vertical lines in the other two panels indicate where x — 2rH. Soo/e are also shown as a 
horizontal line in each panel. 
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Fig. 3. — Evolution of planetary mass and accretion rate without disk dissipation. Figures 
in the left column show accretion rate onto the planet as a function of time, those in the cen- 
tral column show planetary mass as a function of time, and those in the right column show 
accretion rate vs planetary mass. Figures in the top row show the dependence on surface den- 
sity (Soo/(M,r-2) = 10-6(solid line), 10-^ 10-^ lO'^, lO'^ (dot-dashed line)), those in the 
middle row show the dependence on scale height {h/vp — 10~^(solid hne), 10~^-^, 10~^, 10~°-^ 
(dotted line)), and those in the bottom row show the dependence on viscosity coefficient 
(i//(r^fip) = 10-^(solidline),10-^10-^10-^ (dotted line)). The standard values of the 
parameters are Eoo/(M*rp2) = 10"^, i//(r^Op) = 10"^, and h/vp = 10"^ 
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Fig. 4. — Evolution of planetary mass and accretion rate with disk dissipation when h/vp = 
0.1, u/lr^Qp) = 10~^. Three panels are the same axis with Fig. [31 Three thick lines show 
different initial surface densities (Sinit/(M*rp ^) = 10~^'^ (solid line), 10~^'^ (long-dashed line), 
10~^'^ (short-dashed line)) and evolutions without disk dissipation are also shown as thin 
lines. 




Fig. 5. — Final mass of planets as a function of ( (defined at eq.[2l]) in the case with h/vp = 
0.1, i>/{rpQp) = 10~^. The three lines are different initial planetary mass (Mpinit/M* = 
1 X 10-^3 X 10-^l X 10"^). 
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Fig. 6. — Contour lines of log(Mfinai/M*) as a function of C and v/ (rpfip). Left panel is when 
h — O.lrp, right panel is when h — 10~^-^rp. The initial mass is set to 10~^M*. 
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Fig. 7. — Contour lines of log(Mfinai/^*) when global viscous evolution is considered in 
the case with Tdcp = 10^. Horizontal axis is normalized surface density times 10^ (i.e., it 
corresponds to C, when Eoo,init = Sss(rp,0)) and corresponding disk mass assuming R^nt = 
lOrp is also shown on the top. Vertical axis is normalized viscosity and corresponding a is 
also shown on the right. Left panel shows the case with h/r — 0.1 and right panel shows 
when h/r = 0.032. 
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Fig. 8. — Final mass of planets as a function of semi-major axis in AU. Solid line shows 
standard case: a = 0.01, h/r = 0.032 at lAU, r^ep = lO^r, M^isk = 1-3 x lO-^M^, R^^t = 
lOOAU. Dashed line shows Mp^isk (i-e., mass when all viscous- accreting gas is assumed 
to accrete onto planets), and dotted line shows Mp^ss (i-e., mass when the global viscous 
evolution with self-similar solution is not assumed). Initial mass of the planets is set as 
3.2 X 10~^M*, which corresponds to 10 Earth masses in the solar mass system. 
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Fig. 9. — Final mass of planets as a function of semi-major axis in AU for several disk 
models. Thick solid line shows standard described in Figure [H Thick long-dashed 

line shows Mjisk is 10 times smaller case, thick short-dashed line shows viscous a is 10 times 
smaller than the standard case, thick dotted line shows r^ep is 10 times shorter case, thick 
dot-short-dashed line shows the case when scale height is 3 times higher, and thick dot-long- 
dashed line shows the case when gas accretion rate is reduced by a factor of 10. Thin solid 
line and thin dashed line show the masses determined by the viscous and thermals condition 
[a = 0.01), respectively. 
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Table 1. List of notations 



Variables Meaning Definition 



Mp Planetary mass 

Af* Mass of the central star 

-'^disk Disk mass at the initial condition Jq°° 27rrSss('', 0)dr 

^final Final mass of a planet Eq. 1231 

Mflnal,gap Final mass of a planet in the gap-limiting region Eg. 1361 

Affinal diff Final mass of a planet in the diffusion-limiting region Eqs. 1381 and 1391 

Mtrans Planetary mass when Eacc is {l/e)Eoo in the case of 2rH > aim Ea. lB2l 

A^p,init Initial planetary mass 

A^p.crit Maximum planetary mass that does not violate a geometrical limit for oq. 1171 Eq. 1481 

^p, local Final mass if Mp is assumed to be Mp lo^al Eq. l35l 

Mp disk Final mass if Mp is assumed to be M^isk Eg. 1351 

Mp via Planetary mass when viscous condition for gap opening is fulfilled Eq. 1441 

Mp Accretion rate onto a planet Egs. ll6l or 1321 

Mp local Accretion rate onto a planet using cgs. 1141 . 1161 . and 1261 see ^4.3.11 

Mdisk Radial mass flux of the disk due to viscous evolution Eq. 1301 

Soo Unperturbed surface density at planet's orbit Egs.l4l. ll3l .l21 | . l26l 

Sacc Surface density at the accretion band [x = 2rii) Eg. 1161 

Svis Surface density determined by the balance between viscous and gravitational torques Eq. [4] 

Eji Surface density determined by marginally stable state of the Rayleigh condition Eq. 1131 

Saa Surface density of self-similar solution for viscous evolving disk Eq. 1271 

Soo,init Initial unperturbed surface density at planet's orbit Eg. 1211 

rp Semi-major axis of a planet 

Boundary position between the gap-limiting region and the diffusion-limiting region sec i|4.3.3l 

re Boundary position between the diffusion-limiting region and the no-growth region sec i|4.3.3l 

i Position where Evia is (l/e)ScxD Eqs. [3] or [S] 

Xja Position where the Rayleigh condition is marginally fulfilled for Svis Eg. 1101 

T(jiv Time when the accretion rate is the maximum Eg. 1191 

■'"dep Exponential depletion timescale of a disk Eg. 1211 

Tvia Viscous evolution timescale of a disk Eg. 1291 

^lifetime Effective disk lifetime (i.e., shorter one of Tjop and Tvia) Eg. 1341 

f Normalized parameter (Surface density times disk depletion time) Eq. 1241 



